1.1 Chargement des données sur le prix de l’immobilier par EPCI

library(here)
# On situe le dossier dans lequel se trouve nos données
csv_path <- here("data", "donnees_standr.csv")

immo_df <- read.csv2(csv_path)

#Pour visualiser les données dans le doc
datatable(head(immo_df, 10))

Ce fichier est composé des 10 variables suivantes :

• SIREN : code SIREN de l’EPCI
• prix_med : pris médian par EPCI à la vente (au m2 ?)
• perc_log_vac : % logements vacants
• perc_maison : % maisons
• perc_tiny_log : % petits logements (surface < ?)
• dens_pop : densité de population (nb habitants / km2 ?)
• med_niveau_vis : médiane du niveau de vie
• part_log_suroccup : % logements suroccupés
• part_agri_nb_emploi : % agriculteurs
• part_cadre_profintellec_nbemploi : % cadres et professions intellectuelles

La variable SIREN nous servira de “clé” pour joindre ces données statistiques aux données spatiales, la variable prix_med sera la variable que nous chercherons à expliquer (VD), et toutes les autres seront nos variables explicatives (VI).

1.2 Chargement des données géographiques : les EPCI de France métropolitaine

Ces données proviennent de la base ADMIN-EPXRESS-COG de l’IGN, édition 2022. Le format d’entrée est le shapefile mais nous passerons par une conversion au format sf, ce qui nous permet d’utiliser le package mapsf, pour les prévisualiser :

library(here)
library(sf)

shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)

Reading layer `EPCI' from data source 
  `/home/pierson/Travail/projets/formations/2022/LETG_stats_spatiales/GWR_Rzine/data/EPCI.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1242 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 124110 ymin: 6049642 xmax: 1242428 ymax: 7110453
Projected CRS: RGF93_Lambert_93

library(mapsf)
# pour voir les données géographiques
mf_map(x = epci_sf)

# et la table attributaire correspondante
datatable(head(epci_sf, 10))

1.3 Jointure des données géographiques et tabulaires

Les 2 données n’ont pas le même nombre de lignes :

nrow(immo_df)

[1] 1223

nrow(epci_sf)

[1] 1242

On constate que nos deux jeux de données n’ont pas exactement le même nombre de lignes. En effet, le jeu de données immo_df possède moins de lignes que notre objet sf epci_sf. Cela indique simplement que nous n’avons pas l’indication du prix médian de l’immobilier pour tous les EPCI de France métropolitaine. Il peut être intéressant d’identifier et visualiser les EPCI qui n’ont pas de correspondance dans le tableau de données immo_df, pour ce faire on réalise la jointure on conservant toutes les lignes de epci_sf :

# l'option all.x = TRUE permet de garder toutes les lignes de epci_sf,
# même celles qui n'ont pas de correspondance dans immo_df
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN", all.x = TRUE)
nrow(data_immo)

[1] 1242

# on peut filtrer les données de la jointure pour ne voir que les epci n'ayant pas de correspondance dans le tableau immo_df
datatable(data_immo[is.na(data_immo$prix_med),])

mf_map(x = data_immo[is.na(data_immo$prix_med),])

Cependant, notre VD étant prix_med les lignes vides ne nous intéressent pas, nous ne les conserverons pas car elles pourraient poser problème lors de la réalisation de nos analyses. On refait la jointure en ne gardant que les EPCI ayant une correspondance dans le tableau de données :

data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)

[1] 1223

datatable(head(data_immo, 10))

data_immo_sp <- as(data_immo, "Spatial")

2.1 Voisinage

Nous ne développerons pas ici tout ce qu’est et ce qu’implique la définition d’un voisinage. Pour cela, nous vous renvoyons vers les travaux de XXXX ou l’ouvrage XXXXX.

Ici, il est impirtant de savoir qu’un voisinage peut être de 3 types :

• Basé sur la contiguïté,
• Basé sur la distance
• Basé sur la proximité

Lorsque nous travaillons avec des polygones (comme c’est le cas ici), le plus souvent on va se baser sur une matrice de contiguïté. Il faut encore savoir qu’il existe plusieurs types de voisinages basé sur la contiguïté. Dans un cas classique nous utiliserons celui de type queen. Queen est une référence à la reine des échecs. La reine aux échecs peut se déplacer dans toutes les directions, et ici on va considérer les voisins contigus à notre polygone de tous côtés. Il s’oppose au type rook qui fait référence à la tour, les voisins seront donc définis à partir des mouvements de cette pièce sur l’échiquier (dans toutes les directions sauf en diagonale).

Figure 2.8 du manuel INSEE Codifier la structure de voisinage (insee?)

Heureusement R permet assez simplement de définir notre voisinage.

library(spdep)

# Création de la liste des voisins : avec l'option queen = TRUE, 
# sont considérés comme voisins 2 polygones possédant au moins 1 sommet commun
#help(poly2nb)
neighbours_epci <- poly2nb(data_immo, queen = TRUE) 

# Obtention des coordonnées des centroïdes

coord <- st_coordinates(st_centroid(data_immo))

# cette opération se fait aussi avec un shape 
#shape_nbq <- poly2nb(shape, queen = TRUE)
#coord<-coordinates(shape)

Voici la représentation graphique de notre voisinage.

# Faire un graphe de voisinage
plot(neighbours_epci, coord)

Pour comprendre ce que contient neighbours_epci :

# si on prend le 1er élément de neighbours_epci, on voit qu'il a pour voisins les poygones 62, 74 etc.
neighbours_epci[[1]]

[1]   62   74  338 1135 1136 1137 1140

# ce qu'on peut vérifier sur la carte :
neighbours1 <- data_immo[c(1,62,74,338,1135,1136,1137,1140),]
neighbours1$index <- rownames(neighbours1)
mf_map(x = neighbours1)
mf_label(x = neighbours1, var = "index")

Nous précisons qu’un voisinage peut aussi tout à fait se calculer lorsque l’on a pas de polygones mais simplement des coordonnées (des points). Les matrices de distances sont alors souvent plus adaptées. Pour définir le voisinage il faut utiliser les fonctions knearneigh() et knn2nb()

2.2 Création de la matrice de voisinage

Une fois le voisinage défini nous pouvons créer une matrice de voisinage, qui permettra d’attribuer un poids à chaque voisin.

# la fonction nb2listw attribue des poids à chaque voisin
# par ex. si un polygone a 4 voisins, ils auront chacun un poids de 1/4 = 0.25
#help("nb2listw")
neighbours_epci_w <- nb2listw(neighbours_epci)
# les poids sont stockés dans le 3ème élément de neighbours_epci_w
# par ex. si on veut connaître les poids des voisins du 1er élément :
neighbours_epci_w[[3]][1]

[[1]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571

# cet élément a 7 voisins qui ont donc un poids de 1/7 soit ~0.14

3.1 Exploration des variables

Cette première étape très importante permet d’étudier la distribution de nos données et d’identifier d’éventuels individus extrêmes qui pourrait venir perturber les résultats.

Ici par exemple nous serions en droit de nous poser la question de conserver dans notre jeux de données l’entité spatiale avec un prix médian à plus de 10000 qui se détache très clairement de tous les autres EPCI.

Pour visualiser la distribution d’une variable quantitative l’histogramme est une bonne solution. Pour le réaliser nous utilisons dans ce cas le package plotly qui permet l’interactivité de la figure.

# Distribution de la variable dépendante :
library(plotly)
add_histogram(plot_ly(data_immo, x = ~prix_med))

Il est important de réaliser également pour nos VI cet histogramme que nous venons de faire pour notre VD :

# Distribution des variables indépendantes :
a <- add_histogram(plot_ly(data_immo, x = ~log(perc_log_vac), name = "perc_log_vac"))
b <- add_histogram(plot_ly(data_immo, x = ~log(perc_maison), name = "perc_maison"))
c <- add_histogram(plot_ly(data_immo, x = ~log(perc_tiny_log), name = "perc_tiny_log"))
d <- add_histogram(plot_ly(data_immo, x = ~log(dens_pop), name = "dens_pop"))
e <- add_histogram(plot_ly(data_immo, x = ~log(med_niveau_vis), name = "med_niveau_vis"))
f <- add_histogram(plot_ly(data_immo, x = ~log(part_log_suroccup), name = "part_log_suroccup"))
g <- add_histogram(plot_ly(data_immo, x = ~log(part_agri_nb_emploi), name = "part_agri_nb_emploi"))
h <- add_histogram(plot_ly(data_immo, x = ~log(part_cadre_profintellec_nbemploi), name = "part_cadre_profintellec_nbemploi"))
fig = subplot(a, b, c, d, e, f, g, h, nrows = 2)
fig

3.2 Etude des corrélations

Très rapidement, la corrélation permet d’étudier le lien, la relation entre deux variables. La corrélation repose sur la covariance entre les variables.

Quand 2 variables covarient, un écart à la moyenne d’une variable est accompagné par un écart dans le même sens (corrélation positive) ou dans le sens opposé de l’autre (corrélation négative) pour le même individu de l’autre variable. Dit autrement, lorsque deux variables covarient, pour chaque valeur qui s’écarte de la moyenne, on s’attend à trouver un écart à la moyenne pour l’autre variable.

La conception d’un modèle statistique doit absolument être le fruit d’une réflexion portant sur le choix des variables indépendantes (explicatives) et le choix de la méthode de régression. Et pour définir un modèle de régression certaine règles doivent être respectées.

L’étude des corrélations peut donc apporter une aide précieuse dans cette réflexion. Elle pourra nous aider dans le choix des variables à intégrer au modèle mais dans le même temps de vérifier certaines des conditions de réalisation de notre régression.

Ainsi, une analyse des corrélation pourra vérifier :

• l’existence d’un lien entre les variables indépendantes et la variable à étudier. En effet dans une régression linéaire, il est nécessaire d’avoir une relation linéaire entre la VD et les différentes VI.

• la multicolinéarité des variables indépendantes. Les corrélations ne doivent pas être trop fortes entre nos VI. Un coefficient > 0.7 en valeurs absolues doit entraîner la suppression des variables concernées. Cela peut aussi être vérifié très efficacement avec le VIF (Variance Inflation Factor) mais peut se faire seulement après avoir lancé notre modèle.

• L’absence de corrélation entre les variables explicatives du modèle et les variables externes. En effet, les variables d’influence doivent être incluses dans le modèle (sauf dans le cas où cela induirait une trop grande multicolinéarité).

Pour calculer une matrice de corrélation :

library(correlation)

data_cor <- immo_df %>% select(-SIREN)

immo_cor <- correlation(data_cor, redundant = TRUE)

summary(immo_cor)

# Correlation Matrix (pearson-method)

Parameter                        | part_cadre_profintellec_nbemploi | part_agri_nb_emploi | part_log_suroccup | med_niveau_vis | dens_pop | perc_tiny_log | perc_maison | perc_log_vac | prix_med
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
prix_med                         |                          0.55*** |            -0.45*** |           0.60*** |        0.59*** |  0.52*** |       0.49*** |    -0.60*** |     -0.68*** |  1.00***
perc_log_vac                     |                         -0.31*** |             0.39*** |          -0.28*** |       -0.46*** | -0.21*** |      -0.19*** |     0.35*** |      1.00*** |         
perc_maison                      |                         -0.58*** |             0.54*** |          -0.79*** |       -0.24*** | -0.46*** |      -0.75*** |     1.00*** |              |         
perc_tiny_log                    |                          0.75*** |            -0.51*** |           0.87*** |        0.22*** |  0.62*** |       1.00*** |             |              |         
dens_pop                         |                          0.59*** |            -0.29*** |           0.62*** |        0.18*** |  1.00*** |               |             |              |         
med_niveau_vis                   |                          0.43*** |            -0.33*** |           0.15*** |        1.00*** |          |               |             |              |         
part_log_suroccup                |                          0.65*** |            -0.43*** |           1.00*** |                |          |               |             |              |         
part_agri_nb_emploi              |                         -0.57*** |             1.00*** |                   |                |          |               |             |              |         
part_cadre_profintellec_nbemploi |                          1.00*** |                     |                   |                |          |               |             |              |         

p-value adjustment method: Holm (1979)

# On peut aussi visualiser les corrélations à l'aide d'un corrélogramme

# 1- Méthode simple
# immo_cor %>%
#  summary(redundant = FALSE) %>%
#  plot(type="tile", show_labels =TRUE, show_p = TRUE, digits = 1, size_text=3)

# 2- Méthode plus complexe mais meilleure visualisation
mat_cor_comp <- summary(immo_cor, redundant = TRUE)
# Nom des lignes = valeurs de la première colonne ("Parameter")
rownames(mat_cor_comp ) <- mat_cor_comp[,1]
# Transformation du data.frame en objet matrice (+ suppression première colonne)
mat_cor<- as.matrix(mat_cor_comp[,-1])
# Calcul du nombre total d'individus
nb <- nrow(data_cor)
# Calcul des matrices de p-values et des intervalles de confiance
p.value <- cor_to_p(mat_cor, n = nb, method = "auto")
# Extraction de la matrice des p-value uniquement
p_mat <- p.value$p

library(corrplot)
library(RColorBrewer)

corrplot(mat_cor, 
         p.mat = p_mat, 
         type = "upper", 
         order = "hclust", 
         addCoef.col = "white", 
         tl.col = "gray",
         number.cex = 0.5,
         tl.cex= 1,
         tl.srt = 45, 
         col=brewer.pal(n = 8, name = "PRGn"), 
         sig.level = 0.000001, 
         insig = "blank", 
         diag = FALSE, )

L’études des corrélations nous permet ici de confirmer une relation entre notre variable à expliquer et toutes les variables explicatives définies. Elle met aussi à jour de très fortes multicolinéarités, ce qui va nous poser problème. Dans le cadre de la définition d’un modèle linéaire classique, il faudrait sortir du modèle les variables explicatives qui corrèlent trop fortement. Dans le cadre de cette fiche nous faisons le choix de conserver toutes nos VI, cela fera sens au moment de la GWR.

3.3 Régression linéaire ou Méthode des moindre carrés ordinaire (MCO).

La régression linéaire est un des modèles les plus utilisés en SHS, elle peut être simple (une seule variable explicative) ou multiple (plusieurs VI). Le principe n’est en réalité pas si complexe. La régression linéaire consiste à modéliser la covariation entre une variable à expliquer et une ou des variables explicatives.

Pour ce faire le modèle de régression va chercher à estimer les termes de l’équation de la droite de régression entre la variable à expliquer et la variable explicative. Cette équation va prendre la forme \(f(x)= ax + b + e_i\). Équation qui ressemble à la fonction affine \(f(x)= ax + b\), qui est étudiée en général au collège.

Si l’on cherche à modéliser la covariation entre le prix médian et le pourcentage de logements vacants, l’équation de notre droite de régression ressemblerait à : \(prixmed_i = a*perclogvac_i + b + e_i\).

Où \(a\) est le coefficient associé à la variable du pourcentage de logements vacants, \(b\) est la constante qui apparaîtra sous le nom d’intercept dans les résultats et enfin \(e\) qui lui correspond aux résidus. Les résidus étant ce qui incarne l’écart au modèle.

Traditionnellement on va faire usage d’une régression linéaire lorsque l’on veut prédire les valeurs de notre VD dans les cas où elle n’aurait pas été mesurée. Ou si l’on souhaite comprendre les relations statistiques entre les variables.

Pour lancer notre modèle de régression linéaire avec toutes nos VI voici les lignes de commandes :

# Dans le fonctionnement sur R il est important de stocker la régression dans un objet.
# Pour lancer la régression on va utiliser la fonction lm() dont les 2 lettres sont l'acronyme pour linear model
mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi, 
             data = data_immo)

# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
    dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)

Residuals:
    Min      1Q  Median      3Q     Max 
-1566.8  -220.2   -27.2   174.0  3277.3 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
dens_pop                          173.942     15.272  11.389  < 2e-16 ***
med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 391.8 on 1214 degrees of freedom
Multiple R-squared:  0.7784,    Adjusted R-squared:  0.7769 
F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16

#Pour visualiser les résultats de manière plus agréable on peut aussi utiliser la fonction tbl_regression du package gtsummary.
library(gtsummary)
mod.lm %>%
  tbl_regression(intercept = TRUE)

# On peut également visualiser graphiquement les coefficients des variables explicatives

GGally::ggcoef_model(mod.lm)

3.4 Diagnostic de notre modèle linéaire

library(car)

vif(mod.lm)

                    perc_log_vac                      perc_maison 
                        1.577984                         3.204058 
                   perc_tiny_log                         dens_pop 
                        6.221631                         1.864236 
                  med_niveau_vis                part_log_suroccup 
                        1.532403                         5.977951 
             part_agri_nb_emploi part_cadre_profintellec_nbemploi 
                        1.812419                         3.162576 

# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary

library(gtsummary)

mod.lm %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

library(car)

score_vif <- vif(mod.lm)

# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary

barplot(score_vif, main = "VIF Values", horiz = TRUE, col = "steelblue", las=2)
#ajout du seuil de 4
abline(v = 4, lwd = 3, lty = 2)
# et de la limite de 3
abline(v = 3, lwd = 3, lty = 2)

mod.lm2 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)



summary(mod.lm2)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo)

Residuals:
    Min      1Q  Median      3Q     Max 
-1502.2  -226.3   -42.3   173.5  3535.9 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1592.549     11.597 137.329  < 2e-16 ***
perc_log_vac                     -328.301     13.911 -23.601  < 2e-16 ***
perc_maison                      -100.010     20.507  -4.877 1.22e-06 ***
dens_pop                          161.783     15.750  10.272  < 2e-16 ***
med_niveau_vis                    280.586     14.322  19.591  < 2e-16 ***
part_log_suroccup                 237.069     22.792  10.401  < 2e-16 ***
part_agri_nb_emploi                 2.696     15.370   0.175    0.861    
part_cadre_profintellec_nbemploi  -77.744     19.098  -4.071 4.99e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 405.5 on 1215 degrees of freedom
Multiple R-squared:  0.7623,    Adjusted R-squared:  0.761 
F-statistic: 556.8 on 7 and 1215 DF,  p-value: < 2.2e-16

vif(mod.lm2)

                    perc_log_vac                      perc_maison 
                        1.426783                         3.131469 
                        dens_pop                   med_niveau_vis 
                        1.850756                         1.527378 
               part_log_suroccup              part_agri_nb_emploi 
                        3.874580                         1.762159 
part_cadre_profintellec_nbemploi 
                        2.717634 

library(gtsummary)
mod.lm2 %>%
  tbl_regression(intercept = TRUE) %>% add_vif()

GGally::ggcoef_model(mod.lm2)

# L'argument "both" permeet d'utiliser les deux méthodes : ascendante et ascendante
step(mod.lm2, direction = "both")

Start:  AIC=14696.68
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
    part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi

                                   Df Sum of Sq       RSS   AIC
- part_agri_nb_emploi               1      5060 199817567 14695
<none>                                          199812507 14697
- part_cadre_profintellec_nbemploi  1   2725360 202537867 14711
- perc_maison                       1   3911245 203723752 14718
- dens_pop                          1  17351111 217163618 14796
- part_log_suroccup                 1  17792190 217604698 14799
- med_niveau_vis                    1  63121115 262933622 15030
- perc_log_vac                      1  91601521 291414028 15156

Step:  AIC=14694.71
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + 
    part_log_suroccup + part_cadre_profintellec_nbemploi

                                   Df Sum of Sq       RSS   AIC
<none>                                          199817567 14695
+ part_agri_nb_emploi               1      5060 199812507 14697
- part_cadre_profintellec_nbemploi  1   3211704 203029272 14712
- perc_maison                       1   4155413 203972980 14718
- dens_pop                          1  17567092 217384659 14796
- part_log_suroccup                 1  18007058 217824626 14798
- med_niveau_vis                    1  63117884 262935451 15028
- perc_log_vac                      1  94962190 294779757 15168

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + 
    med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, 
    data = data_immo)

Coefficients:
                     (Intercept)                      perc_log_vac  
                         1592.55                           -327.82  
                     perc_maison                          dens_pop  
                          -99.01                            162.05  
                  med_niveau_vis                 part_log_suroccup  
                          280.55                            237.44  
part_cadre_profintellec_nbemploi  
                          -78.93  

mod.lm3 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, data = data_immo)

res_modlm <- mod.lm$residuals
datatable(as.data.frame(res_modlm))

par(mfrow=c(1,3))
qqPlot(mod.lm) #diagramme quantile-quantile qui permet de vérifier l'ajustement d'une distribution à un modèle théorique, ici loi normal

[1]  36 266

hist(rstudent(mod.lm), breaks = 50, col="darkblue", border="white", main="Analyse visuelle des résidus") # Histogramme pour donner une autre indication sur la normalité
plot(rstudent(mod.lm)) # un graphique pour visualiser l'homoscédasticité des résidus

# Pour étudier la normalité on peut utiliser le test de Shapiro-Wilk
shapiro.test(mod.lm$residuals)

    Shapiro-Wilk normality test

data:  mod.lm$residuals
W = 0.90792, p-value < 2.2e-16

# Pour évaluer l'homoscédasticité on peut utiliser le test de Breusch-Pagan. Le package car propose une fonction pour le réaliser
ncvTest(mod.lm)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 1151.844, Df = 1, p = < 2.22e-16

# Pour visualiser les individus concernés
data_immo[c(36, 266),]

Simple feature collection with 2 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 637302.2 ymin: 6532738 xmax: 1013257 ymax: 6879128
Projected CRS: RGF93_Lambert_93
    CODE_SIREN                       ID                                    NOM
36   200023372 EPCI____0000002150235984 CC de la Vallée de Chamonix-Mont-Blanc
266  200054781 EPCI____0000002150236221               Métropole du Grand Paris
                    NATURE prix_med perc_log_vac perc_maison perc_tiny_log
36  Communauté de communes     6400   -2.3174835   -2.558919   -0.07390879
266              Métropole    10920   -0.7044846   -3.732690    6.43751449
      dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi
36  -0.2437935      0.7430146          1.873292           -1.027161
266 22.7709258      0.9077121          7.345980           -1.068415
    part_cadre_profintellec_nbemploi                       geometry
36                        0.07587777 MULTIPOLYGON (((997329.6 65...
266                       6.30681026 MULTIPOLYGON (((666798.7 68...

# Pour relancer un nouveau modèle sans l'individus le plus extrême. Notez que l'on peut en supprimer plusieurs d'un coup avec subset =-c(36,266)
mod.lmx <- update(mod.lm, subset=-266)

# Etudier le nouveau modèle

summary(mod.lmx)

Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
    dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
    part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)

Residuals:
    Min      1Q  Median      3Q     Max 
-1622.6  -215.0   -37.5   170.0  2767.7 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      1587.246     10.526 150.797  < 2e-16 ***
perc_log_vac                     -302.182     13.317 -22.692  < 2e-16 ***
perc_maison                      -132.159     18.814  -7.024 3.57e-12 ***
perc_tiny_log                    -227.353     26.356  -8.626  < 2e-16 ***
dens_pop                           -4.725     19.978  -0.237    0.813    
med_niveau_vis                    284.947     13.012  21.899  < 2e-16 ***
part_log_suroccup                 404.159     25.701  15.725  < 2e-16 ***
part_agri_nb_emploi               -17.545     14.138  -1.241    0.215    
part_cadre_profintellec_nbemploi   23.384     18.841   1.241    0.215    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 367.8 on 1213 degrees of freedom
Multiple R-squared:  0.7823,    Adjusted R-squared:  0.7809 
F-statistic: 544.9 on 8 and 1213 DF,  p-value: < 2.2e-16

vif(mod.lmx)

                    perc_log_vac                      perc_maison 
                        1.589305                         3.168447 
                   perc_tiny_log                         dens_pop 
                        6.073006                         2.089258 
                  med_niveau_vis                part_log_suroccup 
                        1.531885                         5.726830 
             part_agri_nb_emploi part_cadre_profintellec_nbemploi 
                        1.811307                         3.111740 

par(mfrow=c(1,2))
plot(rstudent(mod.lmx)) 
qqPlot(mod.lmx)

[1]  36 180

# Il est possible de comparer les deux modèles et les coefficients

car::compareCoefs(mod.lm, mod.lmx, pvals = TRUE)

Calls:
1: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
  dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
  part_cadre_profintellec_nbemploi, data = data_immo)
2: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + 
  dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + 
  part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)

                                 Model 1 Model 2
(Intercept)                       1592.9  1587.2
SE                                  11.2    10.5
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
perc_log_vac                      -287.3  -302.2
SE                                  14.1    13.3
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
perc_maison                       -128.3  -132.2
SE                                  20.0    18.8
Pr(>|z|)                         1.6e-10 2.1e-12
                                                
perc_tiny_log                     -261.6  -227.4
SE                                  27.9    26.4
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
dens_pop                          173.94   -4.72
SE                                 15.27   19.98
Pr(>|z|)                         < 2e-16    0.81
                                                
med_niveau_vis                     288.0   284.9
SE                                  13.9    13.0
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
part_log_suroccup                  389.0   404.2
SE                                  27.4    25.7
Pr(>|z|)                         < 2e-16 < 2e-16
                                                
part_agri_nb_emploi                -20.8   -17.5
SE                                  15.1    14.1
Pr(>|z|)                            0.17    0.21
                                                
part_cadre_profintellec_nbemploi   -7.84   23.38
SE                                 19.90   18.84
Pr(>|z|)                            0.69    0.21
                                                

library(spdep)

# test de Moran des résidus de la régression H0: pas d'autocorrélation spatiale

lm.morantest(model = mod.lm, 
             listw = neighbours_epci_w)

    Global Moran I for regression residuals

data:  
model: lm(formula = prix_med ~ perc_log_vac + perc_maison +
perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup +
part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data =
data_immo)
weights: neighbours_epci_w

Moran I statistic standard deviate = 21, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.359707408     -0.003517610      0.000299161 

# Test de Geary H0 pas d'autocorrélation.
#  Attention : Pour avoir le  coefficient il faut faire 1-"Résultat test de Geary" (soit ici le coefficient est 0.67)
# Le coefficient de Geary s'étend de 0 à 2, 1 étant le "0" et signifiant aucune corrélation. Par ailleurs, un score inférieurs à 1 implique une corrélation positive et un score supérieur à 1 une corrélation négative.

geary(x = data_immo$prix_med, 
      listw = neighbours_epci_w,
      n = length(neighbours_epci), 
      n1 = length(neighbours_epci)-1, 
      S0 = Szero(neighbours_epci_w))

$C
[1] 0.3324317

$K
[1] 18.24222

data_immo$res_reg <- mod.lm$residuals

# Définition d'une palette de couleur
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Réalisation de la carte
mf_map(x = data_immo, 
       var = "res_reg", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$res_reg,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Résidus de régression\nlinéaire 'classique'", 
       leg_val_rnd = 1)
mf_title("Résidus modèle lm") #titre

4.1 Niveau global

Pour mesure l’autocorrélation comme nous l’avons vue précédemment les deux outils les plus utilisés sont les test de Moran et de Geary.

L’indice de Moran va considérer les variances et covariances en prenant compte de la différence entre chaque et la moyenne de toutes les observations. L’indice de Geary de son côté prend en compte la différence entre les observations voisines.

Pourquoi l’un plutôt que l’autre ? Il semblerait que Moran soit jugé plus stable et revienne plus souvent dans les articles scientifiques.

Ceci dit le coût de réalisation n’est pas très élevé et rien n’empêche de faire l’un et l’autre pour voir si les résultats sont cohérents entre eux.

Commençons par représenter sur une carte notre variable du prix médian des logements par EPCI :

# La palette de couleur:
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")

# Carte du prix médian 
mf_map(x = data_immo, 
       var = "prix_med", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks=quantile(data_immo$prix_med,seq(0,1, by=1/11)), 
       pal=cols_v1, 
       leg_title = "Prix médian", 
       leg_val_rnd = 1)
mf_title("Prix Médian du logement par EPCI France Métropolitaine") #titre

A l’œil, on voit assez nettement qu’une structure spatiale semble exister.

Vérifions le statistiquement !

library(spdep)

# Pour l'occasion on va standardiser notre prix médian. Cela permettra par la suite de le comparer à d'autres variables si d'autres analyses d'autocorrélation spatiale sont réalisées
data_immo$prix_med_z<- (data_immo$prix_med-mean(data_immo$prix_med))/sd(data_immo$prix_med)

# Test de Geary
# Attention à la lecture particulière des résultats de l'indice de Geary
geary.test(data_immo$prix_med_z, neighbours_epci_w, zero.policy = TRUE, randomisation = FALSE) 

    Geary C test under normality

data:  data_immo$prix_med_z 
weights: neighbours_epci_w 

Geary C statistic standard deviate = 35.972, p-value < 2.2e-16
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic       Expectation          Variance 
     0.3324317242      1.0000000000      0.0003444045 

# Test de Moran :

# On indique dans un premier temps la variable que l'on souhaite analyser
# Puis la matrice de voisinage
# L'argument zero.policy=TRUE permet de préciser que l'on souhaite intégrer à l'analyse les entités spatiales qui n'auraient pas de voisins
# L'argument randomisation=FALSE transmet l'instruction à la fonction que nous supposons que la distribution est normale. Dans le cas contraire on devrait partir sur une solution de type Monte-Carlo qui repose sur la randomisation
moran.test(data_immo$prix_med_z, neighbours_epci_w, zero.policy = TRUE, randomisation = FALSE)

    Moran I test under normality

data:  data_immo$prix_med_z  
weights: neighbours_epci_w    

Moran I statistic standard deviate = 41.143, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     0.7154340217     -0.0008183306      0.0003030652 

# Ce test d'autocorrélation se lit exactement comme un test de corrélation classique. On va donc s’intéresser au signe et à la grandeur du coefficient et à la p-value du test.
# Ici on donc une autocorrélation spatiale positive et importante.

Qu’est ce qu’implique l’existence de cette autocorrélation ?

Comme il l’a été mentionné, on pourrait définir l’autocorrélation spatiale comme la corrélation, positive ou négative, d’une variable avec elle-même du fait de la localisation spatiale des observations.

Cette autocorrélation spatiale peut donc être :

• D’une part, le résultat d’un phénomène inobservé (un aléa) ou qu’on ne peut mesurer qui intervient dans l’espace et qui de fait crée une structure spatiale. Il existe différents phénomènes sociaux de la sorte comme par exemple des phénomènes d’interaction ou de diffusion (comme les phénomènes de diffusion technologique). Ces phénomènes peuvent produire de l’autocorrélation spatiale.
• Dautre part, dans la définition du modèle statistique, la mesure de l’autocorrélation spatiale peut être envisagée comme un outil de diagnostic et de détection d’un “mauvais” (du point de vue statistique) modèle (variables non intégrées spatialement corrélées, erreurs sur le choix de l’échelle à laquelle le phénomène spatial est analysé, etc.)

Pour visualiser rapidement la structure spatiale on peut aussi réaliser un diagramme de Moran qui est complémentaire au test statistique.

moran.plot(data_immo$prix_med_z,neighbours_epci_w, labels = TRUE, xlab="prix medians par epci" , ylab="moyenne du prix médians par epci des voisins")

Comment lire ce diagramme ?

Si nos individus sont répartis complètement aléatoirement dans l’espace alors c’est le signe d’une absence d’autocorrélation, la pente de la droite de régression sera donc de 0. Si on contraire la pente est non nulle, comme c’est le cas ici, c’est bien le signe de la présence d’une autocorrélation spatiale.

Un aspect important de ce diagramme est qu’il permet d’ores et déjà de caractériser nos coefficients d’autocorrélation spatiale. Comme ce graphique est centré sur la moyenne qui est de 0, tous les points à droite de l’axe des ordonnées auront une moyenne > 0 et ceux à gauche < 0. Cette réflexion s’applique également à l’axe des abscisses. Par convention on désigne les individus avec une moyenne > 0 par le terme high et ceux avec une moyenne < 0 par le terme low, au sens de supérieur ou inférieur à la moyenne.

Ainsi on peut découper ce diagramme en 4 quadrants. Les quadrants en haut à droite et en-bas à gauche correspondent aux individus ayant une autocorrélation spatiale positive, c’est-à-dire une valeur proche de celle de leurs voisins.

• Pour le quadrant en haut à droite on parle du quadrant High-High composé d’individus ayant une valeur de la variable plus élevée que la moyenne entourés d’autres individus qui leur ressemblent
• Pour le quadrant en bas à gauche on parle du quadrant Low-Low composé d’individus avec un score plus faible que la moyenne avec des voisins avec un score similaire.
• Le quadrant en bas à droite est considéré comme le quadrant High-Low. On y retrouve des individus avec un score plus élevé que la moyenne sur la variable du prix médian mais avec un voisinage qui ne lui ressemble pas : Autocorrélation spatiale négative mais score élevé à la variable.
• Enfin, en haut à gauche on retrouve à l’inverse les individus avec une valeur du prix médian plus faible que la moyenne et un indice d’autocorrélation spatiale négatif. C’est le quadrant Low-High.

4.2 Niveau local

Les indice de Geary et Moran ont une limite assez importante. Ils partent du principe que le processus spatial étudié serait stationnaire, autrement dit global. Cela veux dire que l’effet de la dimension serait le même dans tout notre espace. Ce qui pose sérieusement question et ce d’autant plus que l’emprise géographique augmente. La compréhension et la réalisation d’autocorrélation au niveau local permet d’avancer par la suite vers la GWR.

C’est Luc Anselin qui va développer cette idée et définir le concept d’indicateur d’autocorrélation spatiale locale. Il parlera de LISA (Local Indicators of Spatial Association).

D’après Luc Anselin, le LISA de chaque individu statistique indique l’intensité du regroupement spatial de valeurs similaires autour de cette individu. Dit légèrement autrement un individu avec un LISA élevé va avoir une concentration autour de lui de voisins avec des valeurs similaires( pour nous par exemple un regroupement d’individus avec un prix particulièrement élevé ou à l’inverse particulièrement bas).

Le LISA le plus utilisé est le I de Moran Local.

L’idée est de faire appel ici au LISA et de compléter le niveau global par une approche locale. On va chercher à la fois à détecter des clusters qui correspondent à un regroupement significatif de valeurs identiques dans une zone particulière, et à repérer des zones de non stationnarité, c’est-à-dire qui ne suivent pas le processus global.

Calcul de LISA sur R avec le I de Moran Local

# Pour ce faire on va utiliser le package rgeoda développé également par Luc Anselin pour réaliser sur R les traitements de GeoDa

library(rgeoda)

# Pour utiliser la fonction local_moran proposé par le package rgeoda 2 pré-requis:

# 1- utiliser la fonction queen_weights du package rgeoda pour calculer une matrice de contiguïté de type queen 
queen_w <- queen_weights(data_immo)
# 2- Sortir la variable à étudier dans un vecteur
prix_med_z = data_immo["prix_med_z"]

lisa <- local_moran(queen_w, prix_med_z)

# Pour visualiser les résultats des LISA il faut les stocker dans des objets ou dans des bases de données pour les représenter 
lisa_colors <- lisa_colors(lisa)
lisa_labels <- lisa_labels(lisa)
lisa_clusters <- lisa_clusters(lisa)
lisa_value <- lisa_values(lisa)
lisa_pvalue<- lisa_pvalues(lisa)

Pour illustrer les Moran locaux on peut réaliser une carte :

## Carte de Moran locaux

data_immo$moranlocalvalue<- lisa_values(lisa)

cols_v2 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272",  "#fb6a4a", "#de2d26")


mf_map(x = data_immo, 
       var = "moranlocalvalue", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       pal=cols_v2, 
       leg_title = "Local Moran", 
       leg_val_rnd = 1)
mf_title("Carte des LISA du prix médian du logement")

Si le score du Moran local est > 0 cela implique que l’on a un regroupement de valeurs similaires, plus faibles ou plus élevées que la moyenne. En revanche si le score est < 0 alors on a un regroupement de valeurs dissimilaires, par exemple des valeurs plus faibles entourés de valeur plus fortes (centre de la Corse).

L’étude des p-value associés est également importante car des LISA significatifs (p-value < 0.05) renvoient à des clusters de valeurs (similaires ou dissimilaires) qui seraient plus marqués que ce l’on peut observer si la répartition spatiale du phénomène était aléatoire.

# Carte des p-value des moran locaux
data_immo$moranlocalpvalue<- lisa_pvalues(lisa)

# Pour plus de lisibilité dans la carte on va faire des classes des p-value
data_immo<- data_immo %>%  mutate(lisapvalue_fac = case_when(moranlocalpvalue <= 0.002 ~ "[0.001;0.002[",
                           moranlocalpvalue <= 0.01 ~ "[0.002;0.01[",
                           moranlocalpvalue <= 0.05 ~ "[0.01;0.05[",
                           TRUE ~ "[0.05;0.5]")) %>%
  mutate(lisapvalue_fac = factor(lisapvalue_fac,
                        levels = c("[0.001;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.5]")))


mypal <- mf_get_pal(n = 4, palette = "Reds")

mf_map(x = data_immo, 
       var = "lisapvalue_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal, 
       leg_title = "P-value Local Moran")
mf_title("Carte des significativité des LISA")

Les regroupements que l’on observe vont pouvoir se rapprocher des 4 types de regroupements que nous avions sur le diagramme de Moran.

Cette carte des clusters est représentée selon une convention particulière.

# En utilisant le package mapsf

data_immo$testmoran <- sapply(lisa_clusters, function(x){return(lisa_labels[[x+1]])})

colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")

mf_map(x = data_immo, 
       var = "testmoran", 
       type = "typo", 
       border = "black", 
       lwd = 0.1, 
       pal= colors,
       val_order = c("Not significant","Low-Low","Low-High","High-Low","High-High"),
       leg_title = "Lisa cluster")
mf_title("LISA clusters")

Si on est en présence d’une autocorrélation spatiale au niveau global, les LISA pourront aussi servir d’indicateur d’instabilité locale. Ils indiquent à la fois des clusters locaux qui vont avoir un impact fort sur le processus spatial global (un score d’autocorrélation locale plus marqué que l’autocorrélation globale) ou à l’inverse des zones que se démarquent du processus global (plus faible autocorrélation).

En revanche, s’il n’y a pas d’autocorrélation spatiale au niveau global les LISA vont nous permettre de détecter des zones où des valeurs semblables se regroupent de façon significative. Ils font émerger des structures au niveau local au sein desquelles les liens entre voisins seront particulièrement importants.

5.1 Calcul de la matrice des distance

La première étape est de calculer la distance entre toutes nos observations. Pour ce faire nous utiliserons la fonction gw.dist().

library(GWmodel)

# Le package GWmodel n'est pas compatible avec le format sf il a besoin d'un shape (contrairement à spgwr qui peut travailler avec un format sf)

# Pour construire la matrice de distances entre centroïdes des EPCI :
dm.calib <- gw.dist(dp.locat = coordinates(data_immo_sp))

5.2 Définition de la bande passante

La bande passante est une distance au-delà de laquelle le poids des observations est considéré comme nul. Le calcul de cette distance est très important car la valeur de la bande passante pourra fortement influencer notre modèle. La définition de la bande passante renvoie à quel type de pondération nous souhaitons appliquer. Heureusement la fonction bw.gwr va choisir pour nous le résultat optimal…

Pour ce faire la fonction va se baser sur un critère statistique que l’utilisateur devra définir : le CV (validation croisée) ou le AIC (Critère d’information d’Akaike). Elle reposera aussi sur un type de noyau qu’il faudra également définir : Gaussien, Exponentiel, Bicarré, Tricube ou encore Boxcar. Enfin, il sera également nécessaire de savoir si ce noyau pourra être adaptatif ou fixe.

Voici quelques informations pour guider nos choix :

• Le critère CV a pour objectif de maximiser le pouvoir prédictif du modèle, le critère AIC va chercher un compromis entre le pouvoir prédictif du modèle et son degré de complexité. En général, le critère AIC est privilégié.
• Avec un noyau fixe l’étendue du noyau est déterminée par la distance au point d’intérêt et il est identique en tout point de l’espace. Un noyau fixe est adapté si la répartition des données est homogène dans l’espace, l’unité de la bande passante sera donc une distance. Avec un noyau adaptatif l’étendue du noyau est déterminée par le nombre de voisins. Il est donc plus adapté à une répartition non homogène, l’unité sera alors le nombre de voisins.

Concernant la forme des noyaux :

• Les noyaux gaussiens et exponentiels vont pondérer toutes les observations avec un poids qui tend vers zéro avec la distance au point estimé.
• Les noyaux bisquare et tricube (dont les forme sont très proches) accordent également aux observations un poids décroissant avec la distance, mais par contre ce poids est nul au delà de la distance définie par la bande passante.
• Le noyau Box-Car traite un phénomène continu de façon discontinue.

Sachant que sur la forme du noyau, il est tout à fait possible de comparer deux pondérations et deux modèles de GWR.

# Définition de la bande passante (bandwidth en anglais) :
bw_g <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "gaussian", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17959.04 
Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17884.14 
Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17798.07 
Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17715.24 
Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17622.76 
Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17526.89 
Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 17422.28 
Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 17353.68 
Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 17283.09 
Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 17261.83 
Adaptive bandwidth (number of nearest neighbours): 23 AICc value: 17250.23 
Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 17247.23 
Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79 
Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79 

bw_g

[1] 19

5.3 Estimation du modèle

Une fois la bande passante définie on peut lancer l’estimation de notre modèle de GWR :

mod.gwr_g <- gwr.robust(data = data_immo_sp, 
                   dMat = dm.calib,
                   bw = bw_g,
                   kernel = "gaussian",
                   filtered = FALSE, # un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local. 2 méthodes ont été définies pour gérer cela. 
                                    # Méthode 1 (argument TRUE) on filtre en fonction des individus standardisés. L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
                                    # Méthode 2 (argument FALSE) on diminue le poids des observations aux résidus élevés.
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Si on souhaite comparer deux modèles car nous avons un doute sur les paramètres c’est tout à fait possible. Par exemple ici nous souhaitons comparer deux formes de noyau :

bw_tri <- bw.gwr(data = data_immo_sp, 
              approach = "AICc", 
              kernel = "tricube", 
              adaptive = TRUE, 
              dMat = dm.calib,
              formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)

Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17701.17 
Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17586.94 
Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17422.51 
Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17294.75 
Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17254.1 
Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17279.33 
Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 17259.05 
Adaptive bandwidth (number of nearest neighbours): 112 AICc value: 17257.17 
Adaptive bandwidth (number of nearest neighbours): 138 AICc value: 17254.8 
Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 17253.75 
Adaptive bandwidth (number of nearest neighbours): 117 AICc value: 17255.04 
Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57 
Adaptive bandwidth (number of nearest neighbours): 126 AICc value: 17254.25 
Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57 

mod.gwr_tri <- gwr.robust(data = data_immo_sp, 
                   dMat = dm.calib,
                   bw = bw_tri,
                   kernel = "gaussian",
                   filtered = FALSE,
                   adaptive = TRUE,
                   formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)


Best_gwr <- cbind(
  rbind(bw_g, bw_tri),
  rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
  rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>% 
  `colnames<-`(c("Nb Voisins","R2","AIC")) %>% 
  `rownames<-`(c("GAUSSIAN","TRICUBE"))


Best_gwr

         Nb Voisins        R2      AIC
GAUSSIAN         19 0.9324328 16849.26
TRICUBE         123 0.8577370 17567.05

Le modèle avec une forme qui a été définie au format gaussien explique un meilleur \(R^2\) et le score d’\(AIC\) est plus faible. Ce modèle est donc plus performant et dans notre situation c’est plutôt sur ce modèle qu’il faut partir.

5.4 Interprétation des premiers résultats

Comme pour le modèle linéaire classique, l’objet qui contient notre GWR est composé de plusieurs éléments. Pour obtenir nos résultats il suffit d’appeler l’objet.

# Pour voir les différent élément qui compose notre modèle de GWR
summary(mod.gwr_g)

              Length Class                    Mode
GW.arguments    11   -none-                   list
GW.diagnostic    8   -none-                   list
lm              14   lm                       list
SDF           1223   SpatialPolygonsDataFrame S4  
timings          5   -none-                   list
this.call       13   -none-                   call
Ftests           0   -none-                   list

# Pour accéder aux résultat
mod.gwr_g

   ***********************************************************************
   *                       Package   GWmodel                             *
   ***********************************************************************
   Program starts at: 2022-09-30 16:48:32 
   Call:
   gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel, 
    adaptive = adaptive, p = p, theta = theta, longlat = longlat, 
    dMat = dMat, F123.test = F123.test, cv = T, W.vect = W.vect)

   Dependent (y) variable:  prix_med
   Independent variables:  perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
   Number of data points: 1223
   ***********************************************************************
   *                    Results of Global Regression                     *
   ***********************************************************************

   Call:
    lm(formula = formula, data = data)

   Residuals:
    Min      1Q  Median      3Q     Max 
-1566.8  -220.2   -27.2   174.0  3277.3 

   Coefficients:
                                    Estimate Std. Error t value Pr(>|t|)    
   (Intercept)                      1592.873     11.204 142.171  < 2e-16 ***
   perc_log_vac                     -287.337     14.134 -20.330  < 2e-16 ***
   perc_maison                      -128.255     20.041  -6.400 2.22e-10 ***
   perc_tiny_log                    -261.556     27.934  -9.363  < 2e-16 ***
   dens_pop                          173.942     15.272  11.389  < 2e-16 ***
   med_niveau_vis                    288.017     13.860  20.781  < 2e-16 ***
   part_log_suroccup                 388.982     27.352  14.221  < 2e-16 ***
   part_agri_nb_emploi               -20.785     15.059  -1.380    0.168    
   part_cadre_profintellec_nbemploi   -7.841     19.904  -0.394    0.694    

   ---Significance stars
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
   Residual standard error: 391.8 on 1214 degrees of freedom
   Multiple R-squared: 0.7784
   Adjusted R-squared: 0.7769 
   F-statistic: 532.9 on 8 and 1214 DF,  p-value: < 2.2e-16 
   ***Extra Diagnostic information
   Residual sum of squares: 186354789
   Sigma(hat): 390.6721
   AIC:  18086.13
   AICc:  18086.31
   BIC:  16985.31
   ***********************************************************************
   *          Results of Geographically Weighted Regression              *
   ***********************************************************************

   *********************Model calibration information*********************
   Kernel function: gaussian 
   Adaptive bandwidth: 19 (number of nearest neighbours)
   Regression points: the same locations as observations are used.
   Distance metric: A distance matrix is specified for this model calibration.

   ****************Summary of GWR coefficient estimates:******************
                                           Min.     1st Qu.      Median
   Intercept                         1109.13394  1309.99018  1516.79508
   perc_log_vac                      -919.82407  -224.60028  -163.40316
   perc_maison                       -898.16365  -258.64481  -110.20592
   perc_tiny_log                    -1210.96882  -206.68578   -92.91874
   dens_pop                          -411.20172    72.82448   217.03593
   med_niveau_vis                      81.29015   287.78992   358.32914
   part_log_suroccup                 -543.68600    42.42839   142.48859
   part_agri_nb_emploi               -498.74706   -65.81443   -32.88823
   part_cadre_profintellec_nbemploi  -347.37935   -48.57329     0.58811
                                        3rd Qu.    Max.
   Intercept                         1696.07367 2330.78
   perc_log_vac                      -113.83881  388.64
   perc_maison                         76.01079  896.95
   perc_tiny_log                       32.81544  773.30
   dens_pop                           268.35651  659.84
   med_niveau_vis                     413.36560  853.98
   part_log_suroccup                  247.33650  700.50
   part_agri_nb_emploi                 -1.15602  120.33
   part_cadre_profintellec_nbemploi    49.93194  388.50
   ************************Diagnostic information*************************
   Number of data points: 1223 
   Effective number of parameters (2trace(S) - trace(S'S)): 320.246 
   Effective degrees of freedom (n-2trace(S) + trace(S'S)): 902.754 
   AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 17201.79 
   AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 16849.26 
   BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 17068.01 
   Residual sum of squares: 56808914 
   R-square value:  0.9324328 
   Adjusted R-square value:  0.9084372 

   ***********************************************************************
   Program stops at: 2022-09-30 16:49:07 

Au niveau des résultats, on a d’abord un rappel complet du modèle linéaire classique. Puis vient ensuite les informations concernant notre GWR.

Ce que l’on peut relever immédiatement c’est que le \(R^2 ajusté\) de la GWR est nettement meilleur que celui de la régression linéaire multiple. On passe de \(77\%\) de variance expliqué à \(91\%\) avec la GWR.

La seconde information qui nous intéresse particulièrement sont les coefficient associés à nos VI. On voit immédiatement qu’ils ne sont pas présentés de la même manière que ceux de la régression linéaire. En effet, chaque VI va avoir des coefficients en fonction du minimum, maximum et des quartiles. Cela permet de rendre compte de la stationnarité de l’effet ou non. Dans notre cas on voit qu’il y a bien une variation et même dans certain cas une inversion des signes. Cela laisse supposer une non stationnarité des effets, c’est-à-dire qu’il y aurait un effet local qui ne suivrait pas l’effet global.

Par exemple, pour le pourcentage de logement vacant avec un coefficient global (modèle linéaire) de \(-287\), quand ce pourcentage augmente le prix médian baisse. En simplifiant le pourcentage baisse d’1% le prix médian augmente de \(287€\) (si l’unité du prix médian est l’€). Dans le cas de la densité de population on a un coefficient global positif donc une relation positive. La densité augmente donc le prix médian augmente. Ici, au global la densité augmente de 1 le prix médian augmente de \(173€\).

Les résultats de la GWR illustrent comment les coefficients varient en fonction des unités spatiales, on est bien sur une approche locale. Gardons l’exemple de la densité de population. Dans les lieux où le prix médian est à son minimum le coefficient est de \(-411\) ; On a donc une relation négative. Dans ces espaces la densité augmente de 1 le prix médian baisse de \(411\). Ensuite on va constater une inversion du signe du coefficient. Ainsi dans les EPCI dans le dernier quartile où le prix médian du logement est le plus élevé (par ex. Paris) le coefficient est positif. A son maximum une augmentation d’1 unité entraîne une augmentation du prix de \(659€\). On a donc très clairement un effet de la densité qui ne sera pas du tout le même en fonction du lieu. Ce qui est aussi très intéressant c’est qu’on peut étudier l’intervalle interquartile. Ainsi, toujours pour la densité, ça implique que pour 50% de nos unités spatiales (EPCI entre quartile 1 et 3) une augmentation d’une unité de la densité va impliquer une augmentation du prix médian entre \(72€\) et \(268€\).

La cartographie va être la meilleure manière de représenter les betas (coefficients) et les différents indicateurs fournis avec la GWR, cela nous permet de décrire plus finement et plus précisément les phénomènes observés.

L’ensemble des données est stocké dans le sous objet SDF de notre modèle. Il contient l’ensemble des informations du modèle associé à chaque donnée spatiale.

On peut le convertir en un dataframe pour le visualiser plus facilement. A l’origine il est au format “SpatialPointsDataFrame”.

# Pour visualiser ce fichier dans R
#View(mod.gwr_g$SDF@data)

#Pour voir à quoi il ressemble dans ce document
datatable(mod.gwr_g$SDF@data)

# Pour voir les variables qui le constituent
names(mod.gwr_g$SDF@data)

 [1] "Intercept"                           "perc_log_vac"                       
 [3] "perc_maison"                         "perc_tiny_log"                      
 [5] "dens_pop"                            "med_niveau_vis"                     
 [7] "part_log_suroccup"                   "part_agri_nb_emploi"                
 [9] "part_cadre_profintellec_nbemploi"    "y"                                  
[11] "yhat"                                "residual"                           
[13] "CV_Score"                            "Stud_residual"                      
[15] "Intercept_SE"                        "perc_log_vac_SE"                    
[17] "perc_maison_SE"                      "perc_tiny_log_SE"                   
[19] "dens_pop_SE"                         "med_niveau_vis_SE"                  
[21] "part_log_suroccup_SE"                "part_agri_nb_emploi_SE"             
[23] "part_cadre_profintellec_nbemploi_SE" "Intercept_TV"                       
[25] "perc_log_vac_TV"                     "perc_maison_TV"                     
[27] "perc_tiny_log_TV"                    "dens_pop_TV"                        
[29] "med_niveau_vis_TV"                   "part_log_suroccup_TV"               
[31] "part_agri_nb_emploi_TV"              "part_cadre_profintellec_nbemploi_TV"
[33] "E_weigts"                            "Local_R2"                           

# Intercept : c'est la constante c'est à dire prix médian de référence
# nom de la variable : estimation du coefficient, du beta associé à la VI en chaque point.
# y : les valeurs de la VD
# yhat : valeur de y prédite.
# residual, Stud_residual : résidu et résidu standardisé
# CV_score : score de validation croisée
# _SE : erreur standard de l’estimation du coefficient
# _TV : t-value de l’estimation du coefficient
# E_weight : poids des observations dans la régression robuste
# Local_R2 : R2 au niveau de chaque unité spatiale

res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr

# carto résidus v2
mf_map(x = data_immo, 
       var = "res_gwr", 
       type = "choro", 
       border = "#ebebeb", 
       lwd = 0.1, 
       breaks = quantile(data_immo$res_gwr, seq(0,1, by=1/11)), 
       pal = cols_v2, 
       leg_title = "Résidus GWR", 
       leg_val_rnd = 1)
mf_title("Résidus GWR")

# On ajoute à data_immo les coefficients

data_immo$agri.coef=mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
data_immo$dens_pop.coef=mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef=mod.gwr_g$SDF$med_niveau_vis
data_immo$logvac.coef=mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef=mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef=mod.gwr_g$SDF$part_log_suroccup
data_immo$cadre.coef=mod.gwr_g$SDF$part_cadre_profintellec_nbemploi

# Réaliser la collection des cartes

#par(mfrow = c(2, 4))

mf_map(x = data_immo, var = "agri.coef", type = "choro", pal= "Earth")
mf_title("Coefficients Agriculteurs")

mf_map(x = data_immo, var = "perc_maison.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Maison")

mf_map(x = data_immo, var = "dens_pop.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de dens pop")

mf_map(x = data_immo, var = "med_vie.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Médiane niveau de vie")

mf_map(x = data_immo, var = "logvac.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Logements vacants")

mf_map(x = data_immo, var = "tinylog.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Petits logements")

mf_map(x = data_immo, var = "suroccup.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de logement suroccupés")

mf_map(x = data_immo, var = "cadre.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Cadre")

data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV     

#Définir contrib max
df<- as.data.frame(data_immo)
# On passe les t-values en valeurs absolues pour voir la plus grande contribution dans un sens sens ou dans l'autre
data_immo$contribmax<- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")]


par(mfrow = c(1, 1))

# Carte
mf_map(x = data_immo, var = "contribmax", type = "typo", pal= "Zissou 1")
mf_title("Carte des variables contribuant le plus par epci")

data_immo$r2local=mod.gwr_g$SDF$Local_R2

mf_map(x = data_immo, var = "r2local", type = "choro", pal= "Reds")
mf_title("R2 Locaux")

# Pour rappel si on a plus de 200 individus et le t-value > |1.96| on pourra considérer le coefficient comme significatif au seuil de 0.05 (95% chances que ce ne soit pas dû au hasard)

data_immo$nbsignif_t<- rowSums(abs(df[, c(30:37)]) > 1.96)

mf_map(x = data_immo, var = "nbsignif_t", type = "typo", pal= "Reds")
mf_title("Nombre des Betas significatif par EPCI (t-value)")

# Les p-value ne sont pas fournis dans le modèle de la GWR, on pourrait les calculer à partir de t-value et de l'erreur standard mais le package GWmodel propose une fonction pour les obtenir

pvalue<-gwr.t.adjust(mod.gwr_g)

# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p 
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p

df<- as.data.frame(data_immo)
data_immo$nbsignif_p<- rowSums(df[, c(41:48)] < 0.05)

mf_map(x = data_immo, var = "nbsignif_p", type = "typo", pal= "Reds")
mf_title("Nombre des Betas significatifs par EPCI (p-value)")

# Ici nous représenterons les p-value avec un découpage par classe de significativité et seulement les p-value de 2 VI

par(mfrow = c(1, 2))

# Par exemple les p-value des coefficients de la variable part de l'emploi agriculteur
data_immo<- data_immo %>%  mutate(agri.p_fac = case_when(agri.p<= 0.002 ~ "[0;0.002[",
                           agri.p <= 0.01 ~ "[0.002;0.01[",
                           agri.p <= 0.05 ~ "[0.01;0.05[",
                           agri.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(agri.p_fac = factor(agri.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))


mypal2 <- mf_get_pal(n = 5, palette = "SunsetDark")

mf_map(x = data_immo, 
       var = "agri.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la part d'emploi agriculteurs")

# Pour la densité de population

data_immo<- data_immo %>%  mutate(dens.p_fac = case_when(dens.p<= 0.002 ~ "[0;0.002[",
                           dens.p <= 0.01 ~ "[0.002;0.01[",
                           dens.p <= 0.05 ~ "[0.01;0.05[",
                           dens.p <= 0.1 ~ "[0.05;0.1[",
                           TRUE ~ "[0.1;1]")) %>%
  mutate(dens.p_fac = factor(dens.p_fac,
                        levels = c("[0;0.002[", "[0.002;0.01[",
                                   "[0.01;0.05[",
                                  "[0.05;0.1[", "[0.1;1]")))


mypal2 <- mf_get_pal(n = 5, palette = "SunsetDark")

mf_map(x = data_immo, 
       var = "dens.p_fac", 
       type = "typo", 
       border = "grey3", 
       lwd = 0.1, 
       pal=mypal2, 
       leg_title = "Classe P-value")
mf_title("P-value du coefficient de la densité de population")